package scu.maqiang.meshgeneration;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T31;
import scu.maqiang.mesh.Mesh2T3;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.DirectSSolver;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

import java.util.Arrays;

public class NACA0012 {
    public static void main(String[] args) {
        Mesh2T3 mesh = new Mesh2T3().ring(0.5, 1.0, 80, 20);
        mesh.translate(0.5, 0.0);
        FES2T31 fs = new FES2T31(mesh);
        SRMatrix A = new SRMatrix(fs.GetNdof());
        double[] RHS = new double[fs.GetNdof()];
        fs.assembleStiff(new double[]{1.0}, BVPType.COMMON, A);
        ScalarFunc Ybd1 = (xy, label, param) -> {
            double curve1 =   0.17735 * Math.sqrt(xy[0]) - 0.075597 * xy[0] - 0.212836 * xy[0] * xy[0]
                            + 0.17363 * Math.pow(xy[0], 3) - 0.06254 * Math.pow(xy[0], 4);
            double curve2 = -(0.17735 * Math.sqrt(xy[0]) - 0.075597 * xy[0] - 0.212836 * xy[0] * xy[0]
                            + 0.17363 * Math.pow(xy[0], 3) - 0.06254 * Math.pow(xy[0], 4));
            return xy[1] > 0.0?curve1:curve2;
        };
        ScalarFunc Ybd2 = (xy, label, param) -> xy[1];
        fs.applyBC_MBN(A, Direct.All, 1, 2);
        fs.applyBC_MBN(RHS, Direct.All, new ScalarFunc[]{Ybd1}, null, 1);
        fs.applyBC_MBN(RHS, Direct.All, new ScalarFunc[]{Ybd2}, null, 2);
        DirectSSolver dss = new DirectSSolver(A);
        double[] NACA_Y = dss.CholeskySolve(RHS, 1);
        mesh.setNodeCoordinates(Direct.Y, NACA_Y);
        mesh.toTecplot("NACA0012Mesh.dat");

        A.clear();
        fs.assembleStiff(new double[]{1.0}, BVPType.COMMON, A);
        fs.applyBC_MBN(A, Direct.All, 1, 2);
        fs.applyBC_MBN(RHS, Direct.All, 0.0, 1);
        dss.setMatrix(A);
        NACA_Y = dss.CholeskySolve(RHS, 1);
        double[][] grad = new double[2][fs.GetNdof()];
        fs.computeGradient(NACA_Y, grad);
        double[] ux = MVO.row(grad, 1);
        double[] uy = MVO.row(grad, 0);
        Arrays.setAll(uy, i-> -uy[i]);
        double[] pressure = new double[fs.GetNdof()];
        Arrays.setAll(pressure, i->0.5 * (ux[i] * ux[i] + uy[i] * uy[i]));
        double[][] flowProblem = new double[][]{NACA_Y, ux, uy, pressure};
        mesh.toTecplot("FlowResult.dat", flowProblem);

    }
}
